Band structure, elementary excitations, and stability of a Bose-Einstein condensate in 

a periodic potential 
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We investigate the band structure of a Bose-Einstein condensate in a one-dimensional periodic 
potential by calculating stationary solutions of the Gross-Pitaevskii equation which have the form 
of Bloch waves. We demonstrate that loops ("swallow tails") in the band structure occur both at 
the Brillouin zone boundary and at the center of the zone, and they are therefore a generic feature. 
A physical interpretation of the swallow tails in terms of periodic solitons is given. The linear 
stability of the solutions is investigated as a function of the strength of the mean-field interaction, 
the magnitude of the periodic potential, and the wave vector of the condensate. The regions of 
energetic and dynamical stability are identified by considering the behavior of the Gross-Pitaevskii 
energy functional for small deviations of the condensate wave function from a stationary state. It 
is also shown how for long-wavelength disturbances the stability criteria may be obtained within a 
hydrodynamic approach. 



I. INTRODUCTION 

The possibility of studying experimentally the proper- 
ties of quantum atomic gases in the periodic potential 
created by interference of two laser beams, a so-called 
optical lattice, has led to an explosion of activity, both 
experimental and theoretical. One of the basic proper- 
ties of an atom moving in a periodic potential is that it 
undergoes Bloch oscillations when subjected to a suffi- 
ciently weak external force. In a gas of thermal cesium 
atoms such oscillations were observed as long ago as 1996 
0. Subsequently, a number of experimental studies of 
Bose-Einstein condensates in optical lattices have been 
made. In one-dimensional lattices, interference has been 
observed between condensates initially trapped in differ- 
ent local minima of the potential 0. Also Bloch oscil- 
lations and Josephson oscillations Q of a condensate 
have been observed, and acceleration and collective be- 
havior of a condensate have been studied 5]. In higher- 
dimensional lattices, interference effects have been inves- 
tigated and the transition to an insulating state has been 
observed 

On the theoretical side, Bloch oscillations have been 
investigated Q, and states of uniform flow have been 
studied H, . For sufficiently weak interparticle inter- 
action, the properties of a Bose-Einstein condensate re- 
semble those of a single particle moving in a periodic 
potential. Properties of a condensate in this regime have 
been explored in Ref. [lfj- One of the surprising dis- 
coveries is that the interaction between particles can in- 
fluence the band structure dramatically. In a two-state 
model, Wu and Niu found a loop in the band structure 
at the boundary of the first Brillouin zone |9j and, in 
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detailed calculations of band structure, Wu et al. found 
evidence for non-analytic behavior at the zone boundary 
[Tl| . Another unexpected discovery was that for suffi- 
ciently strong particle interactions, there exists a simple 
exact solution to the Gross-Pitaevskii equation for a con- 
densate with a wave vector k corresponding to the bound- 
ary of the first Brillouin zone 0, 0] . More recently Di- 
akonov et al. have carried out explicit numerical calcula- 
tions of the band structure, and have demonstrated that 
the band has a swallow-tail feature at the zone boundary 
[l4| . They also showed that this behavior is predicted by 
the simple two-component model used previously by Wu 
and Niu @. 

One purpose of this article is to calculate properties 
of stationary states of a Bose-Einstein condensate in an 
optical lattice. We carry out numerical calculations of 
the band structure and investigate the size of the loop 
at the boundary of the first Brillouin zone. In addition, 
it is demonstrated that a similar swallow-tail structure 
can arise at the zone center. It is remarkable that the 
width of the swallow tails remains nonzero even in the 
absence of the lattice potential. We show how this may 
be understood in terms of periodic solitons which, for 
the Gross-Pitaevskii equation, were first investigated by 
Tsuzuki 

A second purpose of the paper is to explore elemen- 
tary excitations and stability of states of uniform super- 
flow. For relatively weak interparticle interactions, this 
has been done in Refs. and ^(j> an d we shall in this 
paper devote most attention to the range of parameters 
for which there are loops in the band structure. Quite 
generally, linear stability of states may be investigated by 
expanding the Gross-Pitaevskii energy functional to sec- 
ond order in the deviation of the condensate wave func- 
tion from the solution for a stationary state. One may 
distinguish two types of stability. The first is energetic 
stability, which is referred to in Ref. 01 as "Landau sta- 
bility" , and the condition for this is that the changes in 
the Gross-Pitaevskii energy functional due to the change 
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in the condensate wave function be positive definite. The 
other form of stability is dynamical stability, and the 
criterion for this is that the linearized time-dependent 
Gross-Pitaevskii equation have no complex eigenvalues. 
To describe excitations with wavelengths long compared 
with the period of the lattice, a hydrodynamic approach 
may be applied. This work represents a generalization 
to moving condensates of the calculations of Ref. 0] 
for condensates at rest. It yields a stability criterion for 
creation of long- wavelength phonons which reduces for a 
translationally-invariant system to the Landau criterion. 

The results obtained in this work are derived from nu- 
merical solutions to the one-dimensional Gross-Pitaevskii 
equation. The method we adopt is to expand the wave 
function in a Fourier series. In order to elucidate the 
physical meaning of our results, we have also carried out 
approximate analytic calculations which yield simple re- 
sults in qualitative agreement with those of the numerical 
calculations. 

In this paper we shall consider only extended states 
having the form of Bloch waves. Due to the nonlinear 
nature of the Gross-Pitaevskii equation there are also 
stationary states corresponding to localized excitations 
such as isolated solitons , as well as states in which the 
density varies periodically with a period different from 
that of the optical potential. 

This paper is organized as follows. In Sec. [H] the prop- 
erties of stationary states are described. There we present 
analytical and numerical calculations, and describe how 
swallow-tail structures may be understood in terms of 
periodic solitons. Section ITTT1 gives a general discussion 
of elementary excitations of a condensate, and energetic 
and dynamical stability. It describes the hydrodynamic 
approach applicable at long wavelengths. Numerical re- 
sults for the stability of a condensate in a one-dimensional 
optical lattice are given in Sec. IIVI Section Ivl contains a 
discussion of our results and concluding remarks. 



II. BLOCH WAVES 

The basic assumption that we shall make in this paper 
is that fluctuation effects are so small that the state of 
the condensate may be calculated in the Gross-Pitaevskii 
approach, in terms of the condensate wave function ip(r). 
The energy of the state is then given by 



dr\ — |V^| 2 + F(r)H 2 
Ira 



1 



(1) 



Here m is the mass of a particle, V(r) is the external 
potential and Uq is the effective interaction between two 
atoms, which is given in terms of the scattering length a 
for two-body collisions by Uo = 4nh 2 a/m. 

Stationary states of the condensate may be found in 
the usual way by demanding that the energy functional 
Q be stationary under variations of ^(r), subject to 



the condition that the total number of particles remain 
unchanged. This yields the time-independent Gross- 
Pitaevskii equation 



2m 



VV + V(r)i> + U \ip\ 2 ip = flip, 



(2) 



where fi is the chemical potential. 

A one-dimensional optical lattice gives rise to a poten- 
tial acting on an atom which has the form 

V{x) = 2V cos 2 (7ra/d) = V cos(2nx/d) + V , (3) 

where d is the period of the lattice. The coefficient Vq, 
which measures the strength of the potential, depends 
on the polarizability of the atom and the intensity of 
the radiation that generates the optical lattice. In fu- 
ture we shall generally neglect the constant term, and 
take the potential to be simply Vq cos(2ttx / d) [25|. In 
addition, we shall not take the potential due to the trap 
into account. Such an approach should give a good ap- 
proximation to the local properties of stationary states of 
the condensate, provided the average density and average 
wave number of the condensate vary slowly in space on 
the scale of the period of the optical lattice. For study- 
ing excitations, this approach will be valid provided the 
wavelength of the excitations is large compared with the 
lattice spacing but small compared with the distance over 
which properties of the unperturbed condensate vary sig- 
nificantly A further assumption we shall make is that the 
states are uniform in the y and z directions. The result- 
ing Gross-Pitaevskii equation has a variety of different 
sorts of stationary solution. Some of these are extended, 
while others, such as solitons, are localized in space. In 
this paper we shall focus on extended solutions to Eq. J2J. 
These are the analogs of Bloch states for a single particle 
in a lattice. As remarked in the introduction, there are 
stationary solutions of the Gross-Pitaevskii equation for 
which the particle density does not have the same period 
as the lattice. As will be explained elsewhere, they are 
related to the self-trapped states of a condensate in a 
double- well potential [H E . Here we shall con- 

fine our attention to solutions of the usual Bloch form 



ip(x) 



Jlx.r 



(4) 



where hk is the quasimomentum and f{x) has the same 
period as the lattice, f(x) = f(x + d). 

The energy per unit volume, E, is then given by 



E 



rd/2 

/ dx 


" h 2 


dip 


/-d/2 


2ra 


dx 



Vq COS 



/ 2ttx 
V~d~ 



ivf + ^wi 4 



(5) 



We determine equilibrium solutions ip of the time- 
independent Gross-Pitaevskii equation by expanding ip 
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in plane waves, 

^max 

V> fe = V^e ite £ a»e i2mx/d , (6) 

' y — — '"'max 

where f is an integer. Here 

i rd/2 

n=~ dx\?P\ 2 (7) 

is the average particle density. From this it follows that 
the coefficients a v satisfy the normalization condition 

'-'max 

£ kl a = i- (8) 

V— — i/ max 

The stationary states of the system may be obtained 
by a variational method, by requiring that the deriva- 
tives of the energy functional JHJ with respect to a„ van- 
ish. There are 2i/ max + 1 complex variables, and one real 
constraint (JSJ. In addition, the overall phase of the wave 
function is arbitrary, so the energy functional depends on 
4^max independent real variables. 

A considerable simplification of the computational ef- 
fort is possible because it turns out that for the station- 
ary states of interest for the range of parameters we have 
considered, the phases <j> v of the coefficients a„ may be 
taken to be the same, or to differ by tt. It is easy to 
demonstrate that such states are indeed stationary un- 
der variations of the phases because the phases occur in 
the energy functional as terms of the type cos(0 vi — <j) V2 ) 
and cos(0 Vl -I- <f> V2 — <f> V3 — <j> Vi )- The derivatives of these 
functions with respect to the phases thus give terms of 
the type sin(0 l/1 - <^„ 2 ) and sin(^ 1 -I- 4> U2 - (f> U3 — 0„ 4 ), 
which vanish if the phases are zero or w. We sought solu- 
tions with other phases by allowing the coefficients to be 
complex, but found none. Because the overall phase is 
arbitrary, this means that we may take the a„ to be real. 
When this is done, there are only 2i/ max independent real 
variables. 

We shall present our results in terms of the energy 
per particle, E/n, as a function of the (one-dimensional) 
wave vector k. As a convenient unit of energy we employ 
the quantity Eq given by 



which is the kinetic energy of a particle with wave vector 
equal to that at the boundary of the first Brillouin zone. 
For an optical lattice made by oppositely directed laser 
beams, the lattice spacing is half the wavelength A of 
the light, and therefore Eq is equal to the kinetic energy 
given to an atom initially at rest when it absorbs a photon 
having the frequency of the laser. 

An unusual feature of the resulting energy bands is the 
appearance of loops in the form of swallow-tail structures, 
as demonstrated in Ref. ^4[. One noteworthy result of 
the present work is that swallow tails can occur also at 



the zone center in higher-lying bands, as illustrated in 
Fig. which shows results of numerical calculations of 
the band structure that will be described in Sec. Ill Bl 
Thus the appearance of swallow tails is a general feature. 
At the zone boundary, the swallow tail appears when the 
interaction energy per particle nUo exceeds the amplitude 
Vo of the potential due to the optical lattice. As the 
parameter nUo grows with respect to Vo, the swallow tail 
increases in width and may extend deep into the zone. 
The condition nUo > Vq is necessary for the swallow-tail 
structure to appear at the zone boundary, k = ir/d. 

Had swallow tails appeared only at the zone bound- 
ary, one might have suspected that their existence was 
related to the fact that for k — ir/d there is an exact 
solution to the Gross-Pitaevskii equation. However, for 
k = there is no exact solution of the Gross-Pitaevskii 
equation, and consequently the existence of swallow tails 
is not connected with the existence of an exact solution. 
As we shall see in detail below, the condition for the ap- 
pearance of a loop at the zone center can be much less 
restrictive than at the zone boundary. Before discussing 
our numerical results further, we shall now analyze a sim- 
ple model for the band structure near the zone center 
which exhibits the main qualitative features of the full 
numerical calculations. 



A. An analytic model 

In Ref. 0| the swallow-tail structure of the lowest en- 
ergy band at the boundary of the first Brillouin zone was 
discussed in terms of a simple trial solution to the Gross- 
Pitaevskii equation. In order to illustrate the generic na- 
ture of the phenomenon we shall here consider the band 
structure at the center of the Brillouin zone, correspond- 
ing to k = 0. 

We employ a trial function of the form 

ip = V^e zkx (a + ai e l2 ™l d + a_ 1 e~ l ^ x ' d ), (10) 

where the coefficients do, a\ and a_i are chosen to be real 
for the reasons given above. The trial function thus mixes 
into the free-particle wave function exp(ikx) states that 
differ by the smallest reciprocal lattice vectors, ±2n/d. 
The normalization condition JHJl is 

al + a\ + a\ = l. (11) 

This constraint is satisfied automatically by expressing 
the coefficients in terms of two angles 9 and <fi according 
to the equations 

ao = cos#, ai=sin#cos0, and a_i = sin 9 sin <j>. 

(12) 

Upon inserting the trial function l|10[l with coefficients 
given by Eq. i|12|) into Eq. (JSJ, we obtain 

E _ 

— — £kin + Cpot + £int, (13) 
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FIG. 1: Energy per particle as a function of wave number for 
the lowest bands. The results are obtained from numerical 
calculations based on the wave function (JSJ, as described in 
SeclUBl 
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FIG. 2: Model calculations of the energy per particle as a 
function of wave number for the same parameters as were 
used in Fig0 The results are obtained using the variational 
method with the trial function 1101 . 



where 



£kir 



2m 



k 2 + 2k ( ^ ) sin 2 6»(cos 2 <f> - sin 2 



2tt 



sin 2 i 



is the kinetic energy, 

e P ot = Vo sin 6 cos #(cos + sin (j)) 
is the potential energy, and 

£int = nUo ( i + cos 2 8 sin 2 9(cos (f> + sin 4>) 

+ ~ sin 4 9 sin 2 2cj) 



(14) 



(15) 



(16) 



is the interaction energy. The stationary points of this en- 
ergy function are obtained by equating to zero the deriva- 
tives of the energy per particle with respect to 6 and <fi. 
The results shown in Fig. [21 were found with this model. 

In order to exhibit the nature of the solutions, we first 
consider the simple case in the absence of interactions 
(U = 0), and with Vo <C E , so the lattice potential 
is a small perturbation. At k — one finds three sta- 
tionary points with different energies. The lowest state 
corresponds to the bottom of the lowest band, and its 
wave function is a plane wave with k = plus a small 
admixture of states with k = ±2n/d. Its energy is given 
to second order in Vo by E = —(V ) 2 /8E . The two 
other states are comprised primarily of plane waves with 
k = ±2ir/d with a small admixture of the k = state. 
The energies of the two states are E = AEq for the state 
which corresponds to the top of the second band, and 



E = AE + (V ) 2 /8E for the state at the bottom of 
the third band. Because of the simplicity of the trial 
function, there are no higher bands in this model. The 
magnitude of the energy gap at k = between the second 
and third bands is thus Vq/8Eq. The reason that it is 
of second order in the lattice potential is that the poten- 
tial is sinusoidal with period d. Consequently, it couples 
directly only states with wave numbers differing by the 
smallest reciprocal lattice vectors, ±2tt/<1 The coupling 
between states with wave vectors k+2n/d and k—2ir/d\s 
indirect, since it is brought about by the coupling of these 
states to the state with wave vector k. In the absence of 
interactions between particles, the stationary points of 
the energy functional for k = have |ai| = |o-i|, that is 
(j> is an odd multiple of 7r/4. 

Figure |3K shows results for non-interacting particles 
(Uq = 0). The familiar band structure is seen. The band 
gap is Vo at k = n/d and V Q 2 /8E at k = for small 
Vo, which is still a good approximation at Vq = 2Eq. 
The lowest band is pushed down in the presence of the 
periodic potential 26]. 

In the presence of interactions (Uo ^ 0), the energy 
landscape in the 6-(j) plane has a more complicated struc- 
ture, and the number of points at which the energy func- 
tional is stationary can be greater than for Uq = Q. Fig- 
ure yp exhibits the resulting band structure. The band 
gap at k = is enhanced by the interaction, in contrast 
to the band gap at k = ir/d, which is reduced. 

Let us now examine the behavior of the band struc- 
ture in the limit of vanishing lattice potential (Vo — * 0). 
As one would expect, the band gaps tend to zero while, 
in contrast, the widths of swallow tails increase with de- 
creasing Vo, and they are nonzero for Vq — > 0. In this 
limit, states corresponding to the upper edge of a swal- 
low tail become degenerate with states in the band above. 
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As we shall discuss in Sec. Ill CI the reason for this sur- 
prising behavior is that states on the upper edge of a 
swallow tail correspond to periodic soliton solutions of 
the Gross-Pitaevskii equation. 

The appearance of a swallow tail requires that the in- 
teraction energy nllo be sufficiently large. In the calcu- 
lations for states near the zone boundary in Ref. 0], it 
was shown that the condition for existence of the swallow 
tail is nllo > Vq. We now carry out a similar calculation 
for the swallow tail at the zone center, using the three- 
state model described above. We investigate the form 
of the energy function (|13H16|1 for k — near the point 
8 = 7r/2, </> = 37r/4. This is a stationary point both for 
Uq = and for [/o / 0. If one imagines that Uo in- 
creases gradually from zero while the other parameters 
remain fixed, one observes in the energy landscape given 
by (|13H16|) that two additional stationary points are split 
off from (9, <fi) — (tt/2, 37r/4) when nllo exceeds a critical 
value. This value may be found by inserting 8 = tt/2 + 5 
and </> = 37r/4 + e into Eqs. ifHfllPoll and expanding the 
energy to second order in 8 and e. The resulting expres- 
sion for the change in the energy per particle relative to 
its value for 5 = e = is 

1-AE - ^nU a ~4E 8 2 + V2V Q Se -^(8 2 + 2e 2 ). 
n 4 2 

(17) 

For nllo = 0, the point 8 = e = is a saddle point. As 
nllo increases, this point turns into a local maximum, and 
two saddle points move out to points with both 6 and e 
nonzero. The condition for the saddle point to turn into 
a local maximum is that the symmetric matrix yielding 
the quadratic form l|17(l have a zero eigenvalue, and this 
occurs for nUo(nUo + 8Eq) = V 2 or 

nU = [(4£ ) 2 + V 2 ] 1/2 - AE . (18) 

For values of the lattice potential small compared with 
Eo, (Vo <C Eq), the critical value obtained from Eq. 118|) 
becomes nllo = V 2 /8Eq. This is physically reasonable, 
since the magnitude of the energy gap in the absence of 
interactions equals V 2 /8Eo, as discussed above. At the 
zone boundary, however, the energy gap in the absence 
of interactions is equal to Vq , corresponding to the con- 
dition nlJo > Vo for the appearance of the swallow-tail 
structure at k = ir/d. Note that the result l|18|) is derived 
from an approximate trial function which becomes exact 
only in the limit Vo <C Eo- Even so, it yields a reason- 
able description of the dependence of the critical value 
of nllo on Vq also for higher values of Vq, as shown by 
the numerical results for the width of the swallow tail in 
Fig. El below. 

The swallow-tail structure illustrated in Figs. 11131 is 
thus a general phenomenon, in that it occurs both at the 
zone boundary and at the zone center, provided the inter- 
action energy nllo is sufficiently large. The phenomenon 
also occurs at other points in the Brillouin zone for solu- 
tions which have a periodicity different from that of the 
optical lattice, as will be demonstrated in future work. 



(a) (b) 
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k k 

FIG. 3: Energy per particle in the first Brillouin zone as 
in Fig [5] The results are obtained by a variational method 
with the trial function given in Eq. 1101 . a) In the absence 
of interaction the band structure (bold curves) exhibits the 
usual band gaps at k — and k = n/d. The band gap is Vo at 
k — Tr/d and Vq /8Eo at k = for small Vo. The thin curves 
show the energies for Vo — » 0, i.e. for the free non-interacting 
system, b) In the presence of interaction the swallow tails 
appear for Uo larger than a critical value, which depends on 
Vo and is different for the two band gaps (bold curves). The 
thin curves illustrate the limit Vo — > 0. 



B. Numerical calculations of band structure 

At the beginning of this section we described the gen- 
eral method used to calculate stationary states of a mov- 
ing condensate. Our numerical procedure is as follows: 
The trial function © is inserted into the energy func- 
tional (jSJ), and stationary points of the resulting expres- 
sion for the energy are found. The normalization of the 
wave function is imposed as in Eq. JSJ, or by generaliz- 
ing Eq. l(T2l to hyperspherical coordinates. We determine 
the solutions by using the Mathematica® routine "Find- 
Root" , which requires as input an initial guess for the 
solution. The latter is found by first solving the problem 
in a reduced basis. Once a solution has been found for a 
particular value of fc, this is used as the initial guess for 
nearby values of k. 

In Fig. ^ two examples of numerical results for the 
band structure are presented. The results are calculated 
with !/ max = 4. This corresponds to a number of basis 
functions equal to 2z/ max + 1. The results for ^ max = 4 
and ^max = 5 differ by less than the thickness of the lines. 
Around k — ir/d the energy of the second band calculated 
with ^ max = 2 differs by less than 1% from the full nu- 
merical result. Around k = and for higher bands, more 
plane waves contribute, cmd typically i^max — 3 is needed 
for 1% convergence. Compared to the simple model used 
in Fig. |21 the numerical results exhibit smaller swallow 
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V Q [units of frVUmd 2 ] 



FIG. 4: Contour plot of the width of the swallow tail in the 
lowest band around k — ir/d. The shaded area indicates the 
region where the swallow tail is absent, i.e., for nUo < Vo- 



tails. The higher-lying energy bands are shifted down- 
wards more than the lower-lying ones, causing the band 
gaps to become narrower. 

Rather than exhibiting the full band structure for dif- 
ferent choices of the parameters, we choose to illustrate 
how one particular feature, the widths of the swallow tails 
at the zone boundary and at the zone center, depends on 
the two dimensionless quantities in the problem, nUo/ Eo 
and Vq/Eq. If the band structure is displayed in the re- 
duced zone scheme, a swallow tail may be split up into 
a number of segments. This is shown in Figs. llTfl where 
the swallow tail at the zone boundary is divided into two 
halves. The full swallow tail may be seen in an extended 
zone representation, and we define the width w of a swal- 
low tail as being the magnitude of the difference between 
the values of k at the two tips of a swallow tail in this 
representation. Widths of swallow tails are exhibited as 
contour plots in Fig. 0] and Fig.|3J 

In Fig. 0] we show a contour plot of the width of the 
swallow tail at the zone boundary in the first band as a 
function of the mean-field interaction nUo and the po- 
tential parameter Vq. The full lines are obtained from 
numerical solutions with ^ max = 3, which is sufficiently 
large to ensure an accuracy better than the thickness of 
the lines. The dotted lines are results for u max = 2. 

The analytic model using the trial function (|10|) over- 
estimates the width by less than 50% for uUq < 8Eo 
and Vq < 8E . Table [I] shows examples of calculated 
widths for different basis sizes. We denote the widths 
calculated for a particular value of i/ max by w Vin ^ . In all 
cases, inclusion of more plane waves improves the result. 
The precision is dependent mainly on the width. Around 
k = Tr/d, fewer plane waves are needed than around k = 



^|"*3 ^1"^: ^1"*= 




V Q [units of frVUmd 2 ] 



FIG. 5: Contour plot of the width of the swallow tail in the 
second band around k — 0. If the condensate wave function 
is approximated by the truncated expression 1101 . this swal- 
low tail is absent if the mean field interaction is less than a 
critical value given by Eq. I|18fl . The corresponding region in 
parameter space is indicated by the shaded area. 



TABLE I: Dependence on basis size of the calculated width 
wv maJI of the swallow tail at the zone boundary. The table 
shows we in units of n/d, and the difference w Uram /w6 — 1 in 
per cent. 









Difference 


Vo nUo 






for f max 








1 


2 


3 


Eo 2E 


0.185547 


13.6 


0.25 


0.00 


E 4E 


0.675781 


27.1 


1.05 


0.02 


Eo 8Eo 


1.56555 


46.7 


3.19 


0.13 


4E 8Eo 


0.552031 


35.0 


3.50 


0.20 


6E 8E 


0.182891 


27.8 


2.69 


0.13 



to describe the wave function to a given accuracy, and 
therefore the precision of the width estimates decreases 
with increasing width for fixed f max - Correspondingly, 
the precision improves with decreasing nUo and increas- 
ing Vq. 

The width of the swallow tail in the second band 
around k — is shown in Fig. As we described in 
the previous subsection, this swallow tail exists provided 
nUo exceeds a critical value, which for small values of 
Vo/Eo is given by Vq /8Eo. The analytic model using the 
trial function I|1U|I estimates the width within 25% for 
nU /E Q < 8 and V /E < 8. Table ITT1 shows some ex- 
amples of the precision within that parameter range. In 
this case, the calculated widths using f max = 1 are closer 
to the full numerical results than are those calculated us- 
ing J/ max = 2. For f max > 2 the precision systematically 
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TABLE II: Dependence on basis size of the calculated width 
ui„ m „ of the swallow tail at the zone center. The table shows 
We in units of n/d, and the difference u>„ max /w6 — 1 in per 
cent. 

Difference 



V nUo 



Eo 4:Eq 
Eo 8Eo 

8Eq &Eq 



0.077969 
0.373789 
0.726953 
0.432695 
0.155234 



for v u 



-0.95 4.36 0.05 0.00 0.00 

3.97 13.4 0.16 0.03 0.00 

10.5 24.0 0.36 0.17 0.01 

0.86 34.0 0.97 0.87 0.06 

-20.3 68.4 3.32 2.52 0.33 



improves. The precision becomes worse as nUo and Vq 
increase. 

The asymptotic behavior of the contours in Fig. for 
Vo — > can be determined analytically within the simple 
model described in Sec. Ill Al The thin lines in Fig. [31 
show an example of the band structure in this limit. For 
fixed nUo, the width of the swallow tail increases with 
decreasing Vo, and it is nonzero for V — 0. At the tip 
of the swallow tail these stationary points of the energy 
function H3J merge at (6,4>) = (vr/2,0) for V -> 0. The 
interaction energy ei„t has a minimum while the kinetic 
energy ekin has a maximum at (7r/2,0). Therefore the 
stationary points merge when the upward curvature of 
eint and the downward curvature of ekin in the (^-direction 
cancel at (n/2, 0): 



d 2 



r £ ki; 



d 2 



rQnt 



6=tt/2 



(19) 



6>=tt/2 



From this criterion we obtain the asymptotic width of 
the swallow tail at the zone center w — (tt / Ad)nUo / Eq 
for Vo — > 0. Compared to the numerical result in Fig. [5] 
this analytic result underestimates the asymptotic value 
of nU for V -> by 7% for w = 0.27r/d, and by 26% for 
w = 2n/d. 

For the swallow tail at the zone boundary, Fig. the 
analysis is further simplified by taking <j> = tt/2, which is 
an exact solution at k = ir/d [LH Il3| . The value 8o at 
which the energy is stationary varies continuously around 
the loop. The tip of the swallow tail is found by setting 
the derivative of k with respect to 9q equal to 0. This 
yields sin2# = —{Va/nUo) 1 ^ 3 , and a width 



TT 

2d 



nUp 
Eo 



2/3 



E 



2/3' 



3/2 



(20) 



Compared to the numerical result in Fig. 0] the contours 
predicted by Eq. H20|l are shifted towards lower values 
of nUo, i.e., the width is somewhat overestimated by 
Eq. JSjjJ. For nU = 2E and Vo = E Q the width given 
by Eq. Q2()|l is 21% above the numerical result and, for 
nU = iE and V = 2E , 38% above. The model is thus 
most precise for small widths. 



C. Physical understanding of swallow tails 

A striking feature of the results above is that the width 
of the swallow tail is nonzero as the strength of the pe- 
riodic potential tends to zero. As we shall now describe, 
the states on the upper edge of the swallow tail then 
correspond to periodic soliton solutions of the Gross- 
Pitaevskii equation first discussed for a condensate in 
the absence of an external potential by Tsuzuki ^{| . For 
potentials having the form of Jacobi elliptic functions, an- 
alytical results for periodic solitons have been obtained 
in Ref. 0. 

The simplest case to think about is that at the zone 
boundary, k = n/d. The solution is then an equally 
spaced array of dark solitons (that is solitons for which 
the wave function vanishes on some surface), with one 
soliton per lattice spacing of the periodic potential. That 
this solution has wave vector ir/d may be seen from the 
fact that the phase difference across a dark soliton is tt. 
Since there is one soliton per lattice spacing, the wave 
vector k, which is the average phase change per unit 
length, is thus ixjd. For states with this value of the 
wave vector, the energy of the highest state in the first 
band and the energy of the state in the second band 
are degenerate if the lattice potential is absent, but the 
degeneracy is broken when a weak periodic potential is 
applied. The state having dark solitons with centers at 
x = rd 7 where r is an integer, has a lower energy, since 
the solitons, which are density rarefactions, are located at 
the maxima of the lattice potential. On the other hand, 
the state having solitons with centers at x = (r + l/2)d 
has a higher energy. 

The situation at k — may be described in similar 
terms. Here the solutions on the upper edge of the swal- 
low tail correspond to periodic solitons with two dark 
solitons for every lattice period. The phase change per 
lattice period is thus 27r, which corresponds to a wave 
vector k = in the reduced zone scheme. When the lat- 
tice potential is applied, there is no change in the energy 
to first order in the lattice potential because neighbor- 
ing dark solitons are separated by d/2 and therefore the 
energy due to the lattice potential vanishes, since the 
potential is purely sinusoidal. However, there will be a 
contribution to second order in Vb since the lattice poten- 
tial will make the spacings between neighboring solitons 
unequal, even-numbered solitons being displaced in one 
direction, and odd-numbered ones in the other direction. 
One may say that the lattice potential causes a dimer- 
ization of the soliton array. This picture gives another 
way of understanding the conclusion arrived at earlier 
that the splitting between a state on the upper edge of 
the swallow tail and that in the next higher band is pro- 
portional to Vo for k = ir/d and proportional to V 2 for 
k = 0. 

Let us now consider states with k ^ 0, ir/d. In the 
absence of a lattice potential, there exist solutions of the 
Gross-Pitaevskii equation which are periodic arrays of 
gray solitons. For these, the density never vanishes, and 
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the phase change across the soliton is less than ir. With 
the boundary conditions that are usually imposed, such 
solitons move with a uniform velocity, which we denote 
by ^soliton- However, by boosting the velocity of the con- 
densate by a constant value — u so iiton everywhere, the so- 
lution becomes a stationary one, and the positions of the 
solitons will remain fixed. The lattice potential will then 
lift the degeneracy of the energy with respect to trans- 
lation of the periodic soliton, just as it did in the ear- 
lier example of an array of dark solitons. The wave vec- 
tor of the condensate is then obtained by combining the 
phase differences due to the solitons with the spatially- 
dependent phase due to the velocity shift, — w S oiiton- As 
the wave vector departs increasingly from k = or ■n / d 1 
the velocity u so iiton and the minimum density in the soli- 
ton increase. Eventually, the density modulation in the 
soliton drops to zero, and the periodic soliton branch 
merges with that for motion of a uniform condensate. If 
the coherence length £ = hj \/2mnUo is much less than 
the lattice spacing, that is nUo ^> E , solitons are gen- 
erally well separated, and the highest value of w S oiiton is 
equal to the sound velocity, (nUo/m) 1 / 2 . 

The physical picture of states corresponding to the 
swallow tail gives insight into the convergence of the wave 
functions as the size of the basis is increased. The charac- 
teristic dimension of an isolated soliton is of order three 
times the coherence length, and therefore the ratio of the 
width of the soliton compared with the lattice spacing 
is of order {Eq/tiUq) 1 / 2 . Thus to give a good account 
of the structure of a soliton, one would expect to need 
i/ max ~ (nUo/Eo) 1 / 2 . For the parameter values that we 
have studied, this is consistent with the fact that the 
numerical calculations were well converged for f max = 3. 



III. ELEMENTARY EXCITATIONS AND 
STABILITY 

In the preceding section we have explored the nature 
of Bloch waves throughout the Brillouin zone, for dif- 
ferent values of the mean-field interaction and the lat- 
tice potential. In the present section we investigate 
their stability to small perturbations. We shall treat 
both energetic and dynamical instabilities, starting from 
the Gross-Pitaevskii energy functional. Subsequently, we 
shall examine the nature of the excitations at long wave- 
lengths, and we shall derive stability criteria from hydro- 
dynamic equations. 



thermodynamic potential G = E — fj,N, where N is the 
particle number, and to allow the variations of ip to be 
arbitrary. Writing tp = tpo + Stp and expanding G to 
second order in Sip, one finds 

G = G[ip ] + SG 1 + SG 2 . (21) 

The first-order term vanishes when ipg satisfies the time- 
independent Gross-Pitaevskii equation J5J. The second- 
order term is 

5G 2 = 

[ dr( 5ip*\-^-V 2 + V(r) - u\Sip (22) 
J \ 2m 

+\u [(r ) 2 (Sip) 2 + ijj 2 (sr) 2 + MM 2 sipsr}^) ■ 

This quadratic form can be written in a compact matrix 
notation as 

SG 2 = ^J drS&ASV, (23) 
where we have introduced the column vector 

»-(£) < 24 » 

and the matrix 

A - In, < 25 » 

\ Wo) L ) 

which is Hermitian. The operator L occurring in (|25|l is 
given by 

L = T + V + 2U \iPo\ 2 -^ (26) 

with T = — H 2 V 2 /2m being the kinetic energy operator. 

The solutions to the Gross-Pitaevskii equation corre- 
spond to stationary values of the thermodynamic poten- 
tial. In order to investigate the stability of these solutions 
under small perturbations we must look at the second or- 
der term (|23l) . When this term is positive for all Sip the 
solution is energetically stable. The system is stable if 
the equation 

A5V = \S^ (27) 

has only positive eigenvalues. Instability sets in when the 
lowest eigenvalue vanishes. 



A. Energetic stability 

To investigate energetic stability of Bloch states, we 
expand the Gross-Pitaevskii energy functional (JIJ to sec- 
ond order in the deviation Sip of the condensate wave 
function from the equilibrium solution ip$, subject to the 
condition that the total particle number be fixed. To 
satisfy the constraint, it is convenient to work with the 



B. Dynamical stability 

To explore dynamical instability, we need to examine 
eigenvalues of the time-dependent Gross-Pitaevskii equa- 
tion, 

h 2 81I1 
-—\7 2 iP + V(r)iP + U \iP\ 2 ilj = ih^ : (28) 
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and its complex conjugate. The pair of equations ob- 
tained by making the substitution ip = ijj + 5tp in this 
equation and its complex conjugate may, when linearized, 
be written in the matrix form 



ih— = a z A5^ 



(29) 




with <x z being the Pauli matrix in the usual representa- 
tion, 



(30) 



Note that a z A is non-Hermitian, and therefore can have 
both real and complex eigenvalues. Complex eigenvalues 
always occur in pairs, since if S^x is an eigenfunction 
of & Z A with eigenvalue A, then a z 5^\ is an eigenfunc- 
tion with eigenvalue A*. Thus if the matrix has com- 
plex eigenvalues, there is always one eigenvalue with a 
positive imaginary part, and therefore the corresponding 
mode will grow exponentially in time. The system is then 
dynamically unstable. 



C. Hydrodynamic analysis 

To study excitations which have wavelengths much 
greater than the lattice spacing, it is possible to use a 
hydrodynamic approach. One works with the average 
particle density n(r) and an average phase to be defined 
below, where the averages are to be taken over a volume 
having linear dimensions much greater than the lattice 
spacing but still much smaller than the wavelength of 
the disturbance. Such an approach has previously been 
employed for small condensate velocities in Ref. [lq. 

Consider a condensate subjected to a potential which 
consists of the sum of two contributions, one due to the 
lattice potential and another one which varies slowly in 
space on the scale of the spacing of the optical lattice. 
We denote the latter contribution by V. If V is spatially 
uniform, the phase </>(r, t) of the condensate wave func- 
tion in a stationary state may be written as the sum of a 
spatially- varying part, <fio(r) and a time-dependent part 

x(t), 



(r,i) =0 o (r)+x(t). 



(31) 



Observe that in writing the phase in this form we have 
nowhere made any assumption about how fast 4>o changes 
over distances of the order of the lattice spacing. The 
phase evolves in time according to the Josephson equa- 
tion 



_dcj> dx 



Of 



dt 



-fi-V, 



(32) 



where /i is the chemical potential calculated for V = 0, 
that is with only the lattice potential acting |27j . 



When the average particle density and the potential V 
vary slowly in space on length scales large compared with 
the lattice spacing, one expects the time rate of change of 
the phase x to t> e given by the same result, except that 
the chemical potential and V now both vary in space. 
In the presence of inhomogeneity, the phase will evolve 
with time at different rates at different points in space, 
thereby "winding up" the phase difference between dif- 
ferent spatial points. The wave vector of the condensate 
wave function is determined by the average rate at which 
the phase of the wave function advances in space. Thus 
the change 6k(r) in the wave vector of the condensate is 
given by 



£k(r) = Vx(r). 



(33) 



It therefore follows from Eq. ilMl'l) that the equation for 
the rate of change of the wave vector has the form 



-V[/i(n,k)+F(r,t)]. 



(34) 



When spatial variations are slow, it is a good approxima- 
tion to assume that the energy density locally is given by 
the expression for the energy density of a uniform sys- 
tem, but with spatially-varying local densities and wave 
vectors: 



E 



dr 



B(n(r),k(r)) +n(r)V 



(35) 



Here E(n, k) is the energy density of the state of the uni- 
form system having a wave vector k and particle density 
n when V = 0. In this approximation, the chemical po- 
tential is that of a bulk system having a wave vector and 
average density equal to the values locally in the nonuni- 
form system. To simplify the notation we omit the bars 
in the following, but it should be remembered that the 
symbols n and k always refer to average values locally. At 
the same level of approximation, the local current den- 
sity, the flux of particle number per unit area, is given 
by the result for a uniform system, 



j = ^V k £(n,k). 
Thus the equation of continuity is 



dn dn 



i 



V • V k £(n,k) = 0. 



(36) 



(37) 



To find the elementary excitations, we now linearize Eqs. 
B34JI and Q37JI ■ We denote changes in the local density 
by Sn, those in the wave vector by <5k, and those in the 
potential by SV. If one looks for solutions varying in 
space and time as expi(q • r — tot), one finds that 



and 



(E nM ■ q - huj)6n + q • E kM ■ 6k = (38) 



E n<n q6n + c\E nM ■ 8k - HtuSk = -qSV. (39) 
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Here 



The derivative 



En n — 



d 2 E 

dn 2 



dn 



E, 



8 2 E 
dkjdki 



nti 2 



(40) 



(41) 



is, apart from factors, a generalization of the usual ef- 
fective mass tensor for a single particle. Note that it 
depends on the particle density and on the wave vector 
of the superfluid flow. The final derivatives are 



En.ki — 



d 2 E 
dndh 



= h 



djj 
dn 



(42) 



In the absence of the lattice potential, it follows from 
Galilean invariance that the contribution to the energy 
per particle that depends on the wave number is given 
by h 2 k 2 /2m. Consequently, the derivative l|4^|l reduces 
to the condensate velocity (times H). All derivatives are 
to be evaluated for the unperturbed value of the density 
n and of the wave vector k. 

The above discussion applies for arbitrary directions 
of k and q. Let us now apply the results to the case 
where both these vectors are in the x direction. The 
eigenfrequencies of the system are found by solving Eqs. 
and with 8V = 0, and they are given by 



hu — q x E n , k ± {E n , n Ek.k.q x 



2\l/2 



(43) 



Equation (|43|) provides a generalization to current- 
carrying states of results derived in Ref. ^(| for a conden- 
sate initially at rest. In order to elucidate the meaning 
of Eq. 1)43(1 . let us first consider the case of k — > 0. The 
energy per particle is quadratic for small fc, and therefore 
E n k tends to zero in this limit and 



fru: = ±(E n . n E kjk q : 



2x1/2 



(44) 



The modes are then sound waves, with the sound speed 
given by 



(45) 



where the derivatives are to be evaluated for k — *■ 0. 
This result agrees with that obtained in Ref. [l6j. For a 
translationally invariant system, E(n, k) — E(n, k = 0) + 
nh 2 k 2 /2m. Therefore = d/.i/dn = Uo and E k , k = 

nH 2 /m, and the sound velocity is thus given by the usual 
result for a homogeneous gas 



nUp 
m 



(46) 



The mixed derivative E n k = h 2 k/m = hv, where v = 
hk/mis the velocity of the fluid. The expression l|43l) for 
the frequency then reduces to the familiar result 22J 



Observe that for a condensate moving in an optical lat- 
tice, the quantity E n ^ k /h = h~ 1 d[i/dk — dj/dn takes 
the place of the mean flow velocity j/n that occurs in 
the analogous result for a translationally invariant sys- 
tem. 

Let us now consider the stability of the system to long- 
wavelength perturbations of the local density and wave 
vector. The system is energetically unstable if such per- 
turbations can lead to a reduction of the energy. In the 
absence of the potential, the functional for the energy 
may be expanded about the original state, and one finds 



E = E + j dr{/i(5n + Hjdk 
E n , n (6n) 2 + 2E, hk SnSk + E k ^8k) 2 ] }. (48) 



The first order terms vanish if the total number of parti- 
cles and the phase of the wave function at the boundaries 
are fixed. The latter condition implies that the change 
in the total particle current is unaltered. The conditions 
for the quadratic form to be positive definite are that 



and 



E n ,n > 0, Ek,k > 



E n ,nE ktk > {E n ,kY 



(49) 



(50) 



Sufficient conditions for energetic stability are that the 
condition (|50|l and one of the conditions (|49|l are satisfied, 
since the other inequality is then satisfied automatically. 
Observe that when condition lj5UI) becomes an equality, 
the system has a zero-frequency mode. 

The numerical calculations to be described in the next 
section indicate that energetic instability sets in first at 
long wavelengths (q — * 0). Consequently, Eqs. I|49|l and 
()50|l are the general conditions for energetic stability. As 
an example, we shall use the condition (| 5 Of) in Sec. IIVI 
below to determine an approximate criterion for the limit 
of stability at the zone boundary. 

The condition for onset of dynamical instability is that 
the eigenfrequency (|43|l become complex, which occurs if 
either E n , n or E k ^ k become negative. The first condition 
corresponds to the compressibility being negative, the 
second to the effective mass being negative. 



IV. STABILITY OF BLOCH WAVES IN A 
ONE-DIMENSIONAL LATTICE 

The stability considerations in the previous section 
were general in nature. In the following we apply them to 
a one-dimensional potential V(x) = Vb cos(27ra;/d) with 
period d. As in Ref. ^l| we consider changes in the con- 
densate wave function of the form 



huj = q • v ± sq. 



(47) 



Si/) 



'ikx 



[u q (x)e^ x +v*Jx)e 



Kj.r I 



(51) 
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where u q (x) and v q (x) have the periodicity of the lattice 
As a result, Eq. I|23[) becomes 

SG 2 = I f drS<S> f B5<S> 



where 



and 



5$ = 



B = 




(52) 



(53) 



(54) 



+ i(±k + q)) +V{x)-^ + 2U \ipo\' i 



with fo(x) = ex-p(—ikx)ipo(x) (cf. Eq. Q}). The opera- 
tors L± is given by 

(d_ 

2m \ dx 

(55) 

According to Eq. I|29[l the linearized time-dependent 
Gross-Pitaevskii equation becomes 



ih- 



dt 



,B5<S>. 



(56) 



As discussed above, the stability of solutions to the 
time-independent Gross-Pitaevskii equation may be de- 
termined by study of the eigenvalues of the operators B 
and a z B. Energetic instability sets in when B first ac- 
quires a zero eigenvalue, while dynamical instability sets 
in when one of the eigenvalues of a z B becomes complex. 
Before presenting numerical results, we give two analyti- 
cal examples of the use of the method. First, we consider 
the case V{x) = 0, and then we derive an approximate 
condition for stability of states at the zone boundary. 



A. The homogeneous Bose gas 



The problem of instability of a homo gene ous Bose gas 
has previously been considered in Ref. [T(|- The calcu- 
lations described in this subsection are similar to those 
of Ref. but we offer a somewhat different physical 
interpretation. 

For a homogeneous gas (V(x) = 0), the solutions to 
the Gross-Pitaevskii equation take the form 

ikx 



'ne 



(57) 



where n is the density. The corresponding chemical po- 
tential is /i = nllo + ti 2 k 2 /2m. The change Sip in the 
condensate wave function is written in the form (|5l)l . 
and since the system is uniform we look for solutions u q 
and v q that are constant in space. The matrix B given 
by Eq. (|54|l then becomes 




2kq) 



(58) 



The stability limit is obtained from the condition that 
the determinant of the matrix vanish, corresponding to 
the existence of a zero eigenvalue. This yields 



(h 2 kq/mf 



(H 2 q 2 /2m ■ 



nU ) 2 - (nU ) 



(59) 



where e q is the Bogoliubov result for the energy of an ex- 
citation in a dilute Bose gas. Thus the condition is equiv- 
alent to the Landau criterion that the minimum velocity 
at which it is energetically favorable to create excitations 
is given by e(q)/hq. On division by q 2 the condition H59|) 
becomes 



7 2 /4- 



(ms/Ky 



(60) 



where the sound velocity s is given by Eq. jloj). 

In addition to the energetic instability considered 
above the system may develop a dynamical instability. 
The dynamical instability exists only when the periodic 
potential is present, since otherwise there is no mecha- 
nism for transferring momentum to the fluid. 

In order to understand the origin of the dynamical in- 
stability for a weak periodic potential let us first consider 
the eigenvalues of a z B in the absence of a periodic po- 
tential. In this case the matrix a z B is obtained from B 
by changing the sign of the matrix elements in the second 
row of Eq. H58|) • Its eigenvalues A are given by 



H 2 kq 



m 



2„2 



± nU 



h z q 



4m 2 / 



1/2 



(61) 



As usual, the physical excitation energies of the system 
correspond to the plus sign in this equation. Note that 
this expression becomes identical with Eq. I|43|) in the 
long- wavelength limit (q — > 0). The eigenvalues l|61|l . 
which are obtained for the case when the periodic poten- 
tial is absent, are always real when Uo is positive. Thus, 
as remarked above, for repulsive interactions there is no 
dynamical instability in the absence of a periodic poten- 
tial. 

Now, let us consider the case of a weak periodic po- 
tential. The appearance of a complex eigenvalue corre- 
sponds to a resonance in which two phonons are created. 
The resonance condition requires the total momentum of 
the phonons to be G = 2n/d (or — G) while their total 
energy must be zero. This implies that 



X+(q) + A+(G - q) = 



(62) 



\k\ = ^ (7 fc oV + 9 4 /4 + \J k 2 (G q) 2 + (G- g)V4 

(63) 

where fco = ms/h. The momenta of the two phonons are 
opposite that of the flow. The condition 163fl is precisely 
the Landau condition for the creation of a pair of excita- 
tions with total wave number G in a superfluid flowing 
with velocity Hk/m. According to Eq. 1)63(1 ■ the mag- 
nitude of the wave number, \k\, for resonance decreases 
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from (Jfeg + G 2 /4) 1/2 to (fc 2 + G 2 /16) 1 / 2 as q increases 
from to G/2. With further increase in q, the magni- 
tude of the resonant wave number increases, since it is 
symmetric with respect to the interchange of q and G — q. 
When the lattice potential is weak, dynamical instability 
therefore first appears when q = G/2 at a wave number 
fc given by |fc| = (fc 2 + G 2 /16) 1 /2. 

In the next section we show how the thresholds for 
energetic and dynamical instability are calculated for a 
non-vanishing periodic potential for different values of 
the wave number fc as functions of nllo and Vq. The 
above prediction for the onset of the dynamical instabil- 
ity for V -> 0, \k\ = {kl + G 2 /!^) 1 ' 2 = (ir/d){nUnj2E + 
1/4) 1 / 2 , is an exact result, first derived in Ref. |lC|. and 
may be compared to the asymptotic value of the numer- 
ical results (solid lines) in Fig. in the limit of Vq — > 0. 
The numerical results agree with the analytical expres- 
sion within the precision of the calculation. 



B. Stability of states at the zone boundary 



obtain the total energy per particle E/n, and for k = n/d 
we get 



„2 y2 

E = jiEo + —Ua — 

2 AU 



and 



d 2 E 
dn 2 



U . 



(69) 



(70) 



The other two derivatives are most easily evaluated by 
using the fact that dE/dk is ft times the particle current 
density (see Eq. (j^El). The current may be calculated 
directly from the trial wave function. For definiteness we 
consider states with positive fc, and the result is 



nh ( 2tt . 2 

fc sin 

a 



(71) 



The derivatives of with respect to n and with respect 
to fc may be calculated by differentiating Eq. |[BSJ, and 
one finds 



A second illustrative example is to consider the con- 
dition for long-wavelength instabilities to arise in a flow 
for which fc = n/d. This may be done using the hydro- 
dynamic formalism described in Sec. IIII CI As shown in 
Refs. 0| and 0, for fc = 7r/rf there is an exact solution 
to the time-independent Gross-Pitaevskii equation of the 
form 



d 2 E ft 2 7T 



lfc *(cos0 + sm0e- i2,n!:/<1 ), 



(64) 



which is the same as Eq. (Uni with <\> = tt/2. To investi- 
gate the stability of the state with fc = n/d, we require 
the solution also for fc in the vicinity of ir/d in order to 
evaluate the derivatives with respect to fc, and we shall 
assume that this is given by Eq. 1)64(1 . with being treated 
as a variational parameter. 

We start from the energy function (|13|) . and for the 
trial function (|64|l the kinetic energy per particle is given 
by 

e ki „ = 4E (k 2 + sin 2 9 - 2k sin 2 6) , (65) 
where k = kd/2-K. The potential energy per particle is 

e P ot = Vosin#cos#, (66) 
while the interaction energy is 
nU 



< U! , = (l + 2 cos 2 9 sin 2 



The energy per particle is stationary when 



(67) 



(k- \) sin 26» = cos 20 + ^ cos 20 sin 20. (68) 

2 oEq oEq 

For fc = n/d, the solution of Eq. j^Hl 1S sin 20 = 
— Vo/nUo- On inserting this result into Eqs. I|65I67|I . we 



dndk md [(nU ) 2 - V 2 } 1 / 2 



and 



8 2 E nh 2 



1 - 



dk 2 m \" nU [(nU ) 2 - V 2 } 



(72) 



(73) 



These derivatives are then inserted into the condition 
lf5T))l. and the boundary of the region of stability is given 
by 



V 2 = (nU ) 



2 nU - 2E Q 
nU + 4£ ' 



(74) 



As we shall see below, this curve deviates by no more 
than 9.5% from the one calculated numerically, i.e., the 
contour for fc = ir/d in Fig. [SJ 



C. Numerical calculations of stability limits 

We now describe results of a stability analysis of the 
stationary state solutions ipo found in Sec. [n] The am- 
plitudes u q and v q in Eq. (|51|l are expanded in terms of 
plane waves: 



i2izlx j d 



l=-l„ 



and 



E 



i2irlx j d 



(75) 



(76) 



l = ~ln 



In the sums, we take ^ max to be less than ^ max . If this is 
not done, spurious instabilities can result. These merely 
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FIG. 6: Contour plot of the maximum wave vector k for ener- 
getic stability. Contours for k < n/d corresponds to states on 
the lowest branch of the lowest band, and those for k > n/d 
to states on the lower edge of the swallow tail at the zone 
boundary. 



express the fact that the condensate wave function has 
not been optimized for modes with wavelengths less than 
d/vmux- The operator B in Eq. (|54H is represented as a 
matrix of dimension 4/ max + 2 in terms of the plane-wave 
basis. 

We now investigate the stability of states correspond- 
ing to points on the swallow tail in the lowest band. In 
the reduced zone representation used in Figs. 11131 this 
swallow tail is split into two or more pieces. However, by 
choosing the first Brillouin zone as < k < 2ir/d, this 
swallow tail will appear in one piece for w < 2n/d. In 
the discussion below, we shall use the latter representa- 
tion, and all values of k will be taken to lie in the interval 
[0, 27r/d]. We shall not consider swallow tails that have a 
width greater than 27t/d. 



1. Energetic stability 

At the boundary of the region of stability, the oper- 
ator B has a vanishing eigenvalue and therefore its de- 
terminant vanishes. We evaluate the conditions under 
which the determinant of B first vanishes by using the 
Mathematica® routine "Det" . The calculation of the en- 
ergetic stability limits proceeds as follows: we choose val- 
ues of k and Vo, and determine the value of uUq at which 
the determinant of B first becomes positive definite for 
all q. We find that energetic instability occurs first for 
q — > 0, in agreement with Fig. 1 of Ref. [Toj . 

FigureEJis a contour plot of the maximum wave vector 
for energetic stability as a function of Vq and uUq. Val- 
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FIG. 7: Contour plot of the maximum wave vector k for 
dynamical stability (solid lines). The notation is the same 
as in Fig. |S] For comparison the corresponding results for 
energetic stability are shown by dotted lines. 

ues of k less than n/d correspond to the lowest energy 
states for a given k, while higher values correspond to 
states on the lower edge of the swallow tail. As the plot 
shows, the range of k values for which states are stable 
increases with increasing uUq and decreasing Vq. The nu- 
merical calculations in Fig. [5] are converged to within the 
thickness of the contour lines for ^ max = 3 and ^ max = 3. 

A number of insights into the behavior of the contours 
may be obtained from analytical arguments. First, the 
intercepts on the wXJq axis of the contours for the wave 
vector at which energetic stability sets in may be deter- 
mined from the Landau criterion. For an interacting Bose 
gas with no lattice, energetic stability occurs when the 
velocity of the gas becomes equal to the sound speed s, 
Eq. l|4^|l . The velocity of the gas is hk/m, and therefore 
the condition is 

h A =s or !^o =2 (^V (77) 

The numerical results agree with this. 

A second general remark is that, for small Vo, one 
would expect on the basis of perturbation theory that 
the contours for energetic stability would behave as Vq, 
again in agreement with the numerical results. However, 
for high values of nUo, the quadratic dependence holds 
only for a limited range of Vq . The contours are approxi- 
mately linear at higher values of Vo, just as are those for 
the width of the swallow tail at the zone boundary (see 
Fig. 0J) . The wave vector of the tip of the swallow tail 
sets a natural limit to the wave vector at which instabil- 
ity sets in, and indeed this limit is approached for large 
nU . 
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A comparison of the width of the swallow tail (Fig. 0} 
and the stability boundary (Fig. EJ) allow us to conclude 
that, within the range of parameters investigated, states 
on the lower edge of the swallow tail at the zone boundary 
are never stable for all k. For the conditions under which 
the spectrum is given by Fig.0i, instability sets in around 
k = l.7n/d and for the conditions appropriate for the 
spectrum shown in Fig.^J), around k = l.lir/d. 

2. Dynamical stability 

The boundary for dynamical stability is determined 
following a numerical procedure similar to the one above 
for energetic stability. On one side of the boundary, the 
operator a z B has only real eigenvalues for all q, while 
on the other side it has some eigenvalues that are com- 
plex. The eigenvalues of a z B are determined using the 
Mathematica® routine "Eigenvalues". With increasing 
k (at fixed nUo and Vb), energetic instability first occurs 
for q — > 0, while dynamical instability first sets in at 
q = ir/d. On the lowest branch of the lowest band, dy- 
namical instability exists only for k > 7r/2<i, in agreement 
with Fig. 1 of Ref. H3- 

Figure \7\ shows a contour plot of the maximum wave 
vector for dynamical stability as a function of Vb and nUo 
as solid lines, and the corresponding contours for ener- 
getic stability are shown as dotted lines. The numerical 
calculations in Fig.0are converged to within the width of 
the contour lines for ^ max = 3 and l max = 3. For given Vq 
and nZ7o, the maximum wave vector for dynamical sta- 
bility is always greater than that for energetic stability. 
The contours of the maximum wave vector for dynamical 
stability are nearly linear for the ranges of Vq and nUo 
investigated. A comparison of the contours for energetic 
and dynamical stability in Fig. shows that the stability 
boundaries almost coincide for larger values of Vq . 



V. DISCUSSION AND CONCLUSIONS 

Our analytical and numerical calculations show that 
the band structure of a Bose condensate in a one- 
dimensional optical lattice is affected dramatically by the 
presence of interactions between particles. The appear- 
ance of swallow-tail-like loop structures is a rather gen- 
eral phenomenon. They occur in the lowest band in the 
vicinity of the zone boundary, as reported earlier |14[ . 
and also, as predicted in this paper, near the zone cen- 
ter in higher bands. Indeed, at band gaps at the zone 
boundary or at the zone center one expects such struc- 
tures to appear quite generally on the band with lower 
energy if the effective interaction between particles is re- 
pulsive, the case studied in this paper. For an attractive 
effective interaction (negative scattering length), the loop 
structures would appear on the upper band at the gap. 
While macroscopic condensates with negative scattering 
length are unstable to collapse if the transverse extent of 



the cloud is large, it might be possible to investigate phe- 
nomena associated with swallow tails in such condensates 
if the condensate is tightly confined in the transverse di- 
rections. 

Analytic results were derived using an approximate 
wave function containing either two plane waves (for 
states with wave vector close to the zone boundary) or 
three plane waves (for states with wave vectors close to 
zero). The coupling between plane wave components in- 
creases with the strengths of the interatomic interaction 
and of the potential Vb, and consequently, the approxi- 
mate wave functions become less accurate. The potential 
induces couplings between plane waves whose wave vec- 
tors differ by only the smallest reciprocal lattice vectors 
±2n/d, while the interaction term introduces couplings 
which are less restricted. The simple analytic expres- 
sions for the energy spectra are in good qualitative and, 
in some cases, quantitative agreement with the numeri- 
cal results. At the zone boundary the approximate wave 
function coincides with the exact one, which accounts for 
the good agreement with the numerical results. 

For swallow tails to appear, the interparticle interac- 
tion must exceed a critical value which depends on the 
band in question. We have derived a simple analytic 
expression, Eq. (|20|) . for the width of the swallow tail 
around the zone boundary as a function of the interac- 
tion nUo and the potential Vb- According to the ana- 
lytical models we have used, in the limit of a vanishing 
potential (Vq — > 0) the width of the swallow tail at the 
zone center behaves as w = nU^ir / AE^d and that at the 
zone boundary as w = tiUqit /2E$d. 

The physical interpretation of states on the upper edge 
of swallow tails is that they are periodic solitons. These 
states exist even in the absence of the lattice potential, 
and this accounts for the fact that the width of the swal- 
low tails does not vanish when the lattice potential is 
absent. 

With respect to experimental observability, an impor- 
tant question is whether or not the stationary states are 
stable. For conditions under which one does not expect 
swallow-tail structures in the lowest band, the stability of 
states has been explored previously by Wu and Niu |10| , 
and for large interactions but only at the zone bound- 
ary in Ref. |ll|. In the present work we have explored 
the stability of states on the swallow tail. States associ- 
ated with the upper edge of the swallow tail are always 
energetically unstable, since they correspond to a saddle 
point in the energy landscape. The energetic and dynam- 
ical stability of the states corresponding to the lower edge 
of the swallow tail have been studied numerically, and we 
found that they become more stable as the strength of 
the interatomic interaction increases. We have not calcu- 
lated growth times for unstable modes, but the calcula- 
tions in Ref. (TJj indicate that these are short compared 
with typical experimental times except under conditions 
very close to the threshold for instability. Thus we ex- 
pect growth of instabilities to be an important effect in 
limiting the conditions under which states corresponding 



15 



to the swallow tail may be investigated experimentally. 

The general formalism for studying stability is rather 
cumbersome, but we showed in Sec. MI Cl that for long- 
wavelength modes one may develop a hydrodynamic ap- 
proach. The detailed numerical calculations indicate that 
energetic stability sets in at long wavelengths and there- 
fore the properties of long-wavelength modes is of direct 
relevance for determining the limit of stability of a con- 
densate. Using the hydrodynamic formalism, we derived 
a simple analytic expression for the stability limit of the 
states at the zone boundary, Eq. 

We have restricted our present study to stationary 
states in which the particle density is periodic with a 
period equal to the spacing d of the lattice. However, as 
will be discussed elsewhere, stationary states with longer 
period, e.g. 2d or 3d, exist. An example of a state 
with a particle density which has a period of two lat- 
tice spacings is a periodic soliton state, with one dark 
soliton for every two lattice cells. The difference in phase 
between two points separated by two lattice spacings, 
4>(x + 2d) — (f>(x) = ±7r, and therefore the wave vector is 
±n/2d. 

In order to observe states corresponding to the lower 
edge of the swallow tail, it is desirable that the states 
be stable. To achieve this would require the mean-field 



energy nUo to be about an order of magnitude larger 
than in current experiments. The criterion for energetic 
stability at k = n/d is nUo > 2Eq, corresponding to a 
a chemical potential /i = 2Eq in the lowest band at the 
zone center, while in the experiment of Cataliotti et al. 
the chemical potential was fi rs 0.2Eq- 

Throughout our calculations we have assumed that 
the system is homogeneous in the directions transverse 
to the optical lattice. However, in actual experiments 
there is usually a confining potential in these directions. 
With sufficiently tight confinement in the transverse di- 
rections, a condensate is expected to behave quasi-one- 
dimensionally, as described in Refs. 0] . Such con- 
densates could provide suitable systems for observing 
some of the non-linear effects predicted in this paper. 
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